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The cluster synchronization (CS) is a very important characteristic for the higher harmonic cou¬ 
pling Kuramoto system. A novel method from the symmetry transformation is provided, and it 
gives CS a profoundly mathematical explanation and clear physical annotation. Detailed numerical 
studies for the order parameters in various conditions confirm the theoretical predictions from this 
new view of the symmetry transformation. The work is very beneficial to the further study on CS 
in various systems. 
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As the simplest and the most celebrated one, the Ku¬ 
ramoto model captures the main property of the collec¬ 
tive synchronization, and is applied in many physical, 
biological and social systems, including electrochemical 
oscillators, Josephson junction arrays, cardiac pacemaker 
cells, circadian rhythms in mammals, network structure 
and neural network[IJ-Q. 

In the large-N limit, the Kuramoto model revealed 
the second continuous transition at the critical coupling 
strength K c . Many generalizations of the Kuramoto 
model have been investigated. Including the large inertia 
in the generalized Kuramoto model, the transition from 
the incoherent state to the collective synchronization be¬ 
came of the first orde®, [f|. Noise also can push the 
incoherent stationary state to become stable 0,0 - When 
the universal coupling strength K becomes oscillator- 
dependent and correlated with the frequency, the explo¬ 
sive synchronization (ES) appears [l]-®]. ES is also an 
abrupt, of the first-like phase transformation. The iden¬ 
tical oscillators with the nonlocal coupling strength will 
give rise to the new chimera state, which is the combi¬ 
nation of the coherent state and the incoherent state for 
the identical oscillators fid)-[l8j] . 


All the above examples are of the first harmonic cou¬ 
pling as H(0j — 0i) = K,j sin(0,j — 0i). Whenever the 
higher globally coupling harmonic term is introduced, in¬ 
teresting phenomena appear, like the cluster synchroniza¬ 
tion (CS) or multi-entrainment, and switching of the os¬ 
cillators between different clusters with the external force 
®-®. The higher harmonic coupling is dominating 


in (f -Josephson junction j21 L 22], in the electrochemical 
oscillators in higher voltage jl9L l23l. l24j , in neuronal net¬ 
works with learning and network adaption®]- [30]. CS 
is the most outstanding feature of the Kuramoto model 
with higher order harmonic coupling. 

CS or multi-entrainment has been investigated by the 
method of self-consistent approach in Refs.®.® ,®- 
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[32| . Neural network actually studied the combination 
of the first and second harmonic couplings in the gener¬ 
alized Kuramoto model®]-®], which is also treated in 
Ref.®, ®. In the N identical oscillators’case, the sym¬ 
metry viewpoint is applied and CS of the two groups of m 
and N — m oscillators is connected with their symmetry 
groups of the dynamics S m x Sw-m 28], ®. The sym¬ 


metry group Sn is only suited for the identical oscillators 
in the Kuramoto model. However, it is still very difficult 
to obtain clear analytical results by the self-consistent 
approach and detailed understanding of the stability of 
the asynchronous states is still missing ®, ®. 

In order to overcome the shortcoming, the den¬ 
sity function for the second harmonic coupling case 
is decomposed into the symmetric and asymmetric 
parts in Ref. , and the Ott-Antonsen (OA) mech¬ 
anism is utilized to analyze the symmetric case. Ott- 
Antonsen ansatz (reduction mechanism) is powerful in 
obtain ing the analytical understanding of the Kuramoto 
model[34| and has been applied in generalized Kuramoto 
ones, like the forced Kuramoto [35], the bimodal fre¬ 
quency distribution®, ®] and the second harmonic 
coupling [19]. The asymmetric clustering is also showed 
in Ref. [191] very sensitive to the non-uniform initial con¬ 
dition. Even the OA ansatz can not give a clear physical 
understanding of CS thanks to unavailability in obtain¬ 
ing the analytic form for the intricate asymmetry density 
function 

Here we will investigate the higher-harmonic coupling 
Kuramoto model from the point of symmetry, and pro¬ 
vide a group transformation, which is completely differ¬ 
ent from that of Sn group. Then we give CS a thoroughly 
novel interpretation. 

In the original Kuramoto model 


1 N 

8 n =u n + — sin (Oj - 6 n ), (1) 


3 =1 


the coupling strength K > 0 is assumed. The higher 
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harmonic coupling of the generalized Kuramoto model is 

1 N 

0 n = w n +— ^2,K rn su\m(6j - 0 n ). ( 2 ) 

JV j =i 

The order parameter is defined as 


re " = vX>"'- P) 

J=1 

Ref. 0 shows that the critical parameter K ^ relating 
with the corresponding m-th order parameter is the same 
K = 2A for all integers m, where 2A is the width of the 
original Lorentz frequency distribution in Ref. 0. In the 
case of small strength K m < K^, the term dominates 
the change of the phase 0 n and the whole phase system is 
in the incoherent state. Whenever K m exceeds K the 
second terms in Eq.® predominate and CS emerges 0 . 

Why the critical coupling K ^ is the same for all dif¬ 
ferent integers ml Is it only a coincidence or is there 
some underlying reason? The study in the letter shows 
that the symmetry that the Kuramoto system keeps is 
responsible for. 

By introduction of the transformation 

<t> = rnO, u} n = mu> n , K = mK m , (4) 

Eq.® takes the form 

1 N 

4n = Un + Jr K sin (<Aj “ 4>n), (5) 

iV J=1 

which is the same as Eq. m of the standard Kuramoto 
model0. 

The generalized order parameters is defined as i?i = 
i e '^ = \Ri\e llfil , and Eq.® becomes </>'„ = 
Q n + — R\e l(t,n ). In the limit N — > oo, the 

density function /(</>, ui,t) is introduced to describe the 
distribution of the phases at a given frequency and sat¬ 
isfies the continuous equations 


d t f + 3^ 





( 6 ) 


The solutions to Eq.® could be obtained from the ones 
to Eq.©, where OA mechanism could be utilized when 
the distribution of the natural frequency is the Lorentz’s 
one |.T1 . 

The density distribution function /(</>, ui,t) is a peri¬ 
odic function satisfying f(<j> + 2n,u),t) = f(cj),u],t), and 
in the case to = 2 it corresponds to the symmetric one 
in Ref0]. It is easy to see that in the stationary state 
f(<f>,u!,t) = /(</>,w,t)|t=oo the system is partially syn¬ 
chronic whenever K > 2A for the Lorentz distribution of 
the natural frequency of the phases with A its width fl9i], 
[13, 0|. Hence the same critical parameter K c is real¬ 
ized for all TO-th order parameters. Combination of the 


symmetry transformation ©, f(<f>,u>,t) can completely 
determine the evolution of the dynamics and this is the 
key point in the latter to study CS. 

We apply both the transformation © and the distribu¬ 
tion function t) with its periodic property to inves¬ 

tigate the order parameters for the oscillators in several 
special cases, and make the corresponding predictions on 
the order parameters. Then we integrate Eq.® directly 
for these cases and obtain the corresponding numerical 
order parameters directly from the numerical integration. 
The later numerical results confirm the former prediction. 
The details are divided into five groups in the following. 

(a)If the initial oscillators’ phases are uniformly dis¬ 
tributed in (0, 2n), together with the transformations ® , 
the order parameter r in Eq. © in the large N limit turns 
out as 

i /*+oo /*2?n7r 

re^ = — / / f((/),uj,t)e l " : du}d(j) = 0, (7) 

mi-oc Jo 

where the upper integral number is 2 to 7 t due to the trans¬ 
formations ©■ The symmetry property of the distribu¬ 
tion function f{<f> + 2n,u),t) = f(<j>,w,t) makes the or¬ 
der parameter r = 0, no matter what a great coupling 
strength is applied. This is actually the manifestation 
of the cluster property of the higher harmonic coupling. 
The TO-term harmonic coupling will give rise to the cor¬ 
responding to clusters, and the phase oscillators in one 
cluster behave completely the same way as those in an¬ 
other cluster. The order parameter r can only take the 
zero value, which is a typical manifestations for the clus¬ 
ter phenomena of the higher harmonic coupling in the 
system. 



FIG. 1: Diagrams r(t ) of initially uniform distribution of the 
phases with the same coupling strength K — 1, K < K c in the 
left panel indicating the incoherent state and K = 5, K > K c 
in the right panel showing CS for m = 2, 3. Notice that the 
green lines are shielded by the red lines. The number of the 
oscillators is N in all the figures 


The numerical studies in FigQ] confirm the above ideas 
of r « 0 with the uniform initial distribution in (0, 27 t) of 
the phases and different coupling strengths and different 
higher-term couplings. When K < K c , all parameters 
r in the three cases in FigJT] is similar, and the oscilla¬ 
tors are incoherent. Of course, the symmetry properties 
in Eq.® will force parameters r in to > 2 will be more 
smaller than that in to = 1 case. When the coupling 
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strengths K surpass K c , r s=a 1 is obtained for usual Ku- 
ramoto model (m=l), and r ss 0 stands out in m > 2 
in Figdl which indicates the formation of the clustering 
synchronization. Fig[5] further gives clustering synchro¬ 
nization for the case of K = 5, m = 2, to = 3 by the 
phases’ position in the circles in the middle and the left 
panels. 



FIG. 2: The cluster synchronization of the phases in the cir¬ 
cles corresponding the parameters in the second panel in Fig.l 
with the same parameters I\ = 5 but different m = 1,2,3. 
There are 1000 oscillators and their positions are indicated 
by the small blue circles in the big circles. 


Furthermore, the more higher harmonic coupling is ap¬ 
plied, the more clusters appear and the smaller section 
of the whole range 2-7T every cluster occupies. Therefore, 
very higher harmonic coupling of the oscillators will pro¬ 
duce psudo-synchronization where the large number of 
the clusters will globally behave the same way with the 
oscillators within each very small cluster being any state. 


In each cluster section, the phase oscillators may stay 
in the incoherent state or partially synchronic state, 
even in synchronization state depending on whether or 
not K > K c . The order parameter for the phase 
oscillators in each cluster is the same and is rf = 

^ 0 for the first clus " 
ter, where </>/ is the center of the first cluster. The order 
parameters rf and r are different from each other, but 
the critical coupling strengthes for them are the same. 
The cause is that both rf and r are determined by the 
same distribution function f(<p,U},t ), which decides the 
critical strength K c . 

(b) When the initial phases distribution (0, A) is nar¬ 
rowed and much smaller than —, and the coupling 
strength is stronger than the critical K c , and with ev¬ 
ery term in the second part ]CyLi K m sin m(6j — 6 n ) 
in Eq.@ has the same effect as the corresponding term 
in ordinary Kuramotto dT]) , and they together dominate 
over the first part and attract all oscillators to the syn¬ 
chronic state. So, the initial synchronized state in one 
cluster will remain synchronized all the time, just like in 
the ordinary Kuramoto model and the order parameter 

r = I f-oo So* /(0> w, « rf « 1 is realized. 

In this case, almost all oscillators are synchronized in 
the only one cluster, and no other cluster synchroniza¬ 
tion exists. See the phases in the circles in the second 
(A = 0.757T, to = 2) and third ( A = 0.87T, to = 2) panels 
in Fig[3]for details. Note there are several phases are left 
opposite to the synchronized one in the two panels. 


J 



FIG. 3: Schematic diagrams r(t) of initial uniform phases distribution in (0, A) with the same parameters K = 5, m = 2 but 
different A. 


(c)When A of the initial distribution (0, A) goes beyond 
2 ^”»i hut less than 2-7T, there approximately emerges m 
cluster synchronization. As the order parameter r « 0 is 
achieved, which is shown in Figl3]with A = 1.2n, 1.257T 
in the case of to = 2, K = 5 for the Gauss frequency 
distribution with K c < 2. In this case, the initial range of 
(f> = m6 exceeds 2 (to — l)7r, and guarantee the validation 
of the transformation © and Eq. ©■ So CS is again 
connected with the symmetry of the Kuramoto model 
©. The fourth (A = 1.257T, to = 2) and fifth (A = 
1.27r, to = 2) panels in Fig[3] agree with this analysis. 


(d)When A of the initial distribution (0, A) lies in 
2(n-i)7r 2 raw w j^j 1 n < m an d rn greater than 2, the 

mm 07 

partial CS appears indicated by the order parameter r 
approaches neither 0 nor 1. The order parameter takes 
the form 

/ +oo nmA 

/ /(</>, w,t)e J ^dujd(/>|, (8) 

-oo J 0 

Generally, in the stationary state, the upper bound toA 
could be replaced by 2mr in most cases. Roughly, 
r(A) ss eZ_ ™~) r ilj where rf is the first clus- 
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ter’s order parameter, and is near 1 in the case K be¬ 
ing much large than K c . The numerical studies in 
Fig© crudely illustrate the feature for r(A). Define 
r(n) = I f-x fo’ l7r /(</>, it is easy to see 

that r(n) < r(A) < r(n — 1). For example, r(n) 
takes the following numerical values r(l) = 1, r( 2) ~ 
0.8666, r( 3) « 0.6777, r(4) « 0.4333, r(5) = 0.2, r(6) = 
0 for the case m = 6 and K > K c . Fig© shows the data 
for r(A) meets the constrain as above mentioned. This 
again tells one that the distribution function /(</>, u>,t) 
can be used to obtain the parameter r(A) through the 
symmetry transformation ©• 



FIG. 4: Detailed diagrams r(t) on the initial uniform phases 
distribution in (0, A) with the same parameters K = 5, m = 6 
but different A. 


(e)When A of the initial distribution (0,A) is , 
the dynamics of the model is very sensitive to the initial 
conditions. The order parameters is defined as 

-i /*+oo /'2nn 

r(A) = -\ / + A|. (9) 

n J -oo Jo 

The oscillators with initial phase very near can either 
lag into the n-th cluster or drift forward to the (n + 
l)-th cluster. A is the related order parameter for the 
entering the (n + l)-th cluster, and the fraction of the 
oscillators for this kind is very sensive to the system’s 
initial conditions. So do the parameters A and r(A). n 
in Eq.© is in (n,n + 1) depending on parameter A. All 
are shown in Fig.??. 

We now conclude the symmetric viewpoint. The trans¬ 
formations © make it possible to relate the results of 
Eq.© with that of Eq.©. Furthermore, the periods for 
both mO in Eq.© and 4> i n Eq.© are identical, that 
is, 27 t. The case K m > K c > 0 in the model © is 
the same case K > K c for the transformed model ©• 
So the initial phases distribution in (0,A) in Eq.© is 
equivalent to the initial distribution (0 ,mA) in Eq.©. 
So when A £ (^L, 2 ^ Tt ^ 1 ^ 7r ), the initial phases distri¬ 
bution for variable </> is (0,2(n + l)7r). Every 27r dis¬ 
tribution in Eq.© will be a Kuramoto model and will 
be synchronized when K > K c . Therefore the distri¬ 
bution (0, 2(n + l)7r) for (f> will be equivalent to n + 1 
Kuramoto models, that is, the n + 1 cluster synchro¬ 
nization. This is the root of the formula ©• Mathe¬ 
matically, the above results means that the phases lines 



FIG. 5: Schematic diagrams r(f) of initial uniform phases 
distribution in (0, ir) with the same parameters K = 5, m = 
2. The only difference for the two cases is their initial values 
for the phases. 

n = 0,1, • • ■ , in — 1 in the unit circle divide the 
interval and the switching across them is forbidden gen¬ 
erally when K > K c . Only the phases very very close 
to them could separate into either forward interval or 
backward interval due to their different frequencies. The 
forbidden lines give vivid demonstration of CS. 

In addition, Eq.© is also invariant under the transla¬ 
tion 9 n = 9 n + 9 q, so CS phenomena is invariant under 
the translation of the initial conditions, see Fig.?? for 
details. 



FIG. 6: Schematic demonstration of the translation properties 
for r(t) of initial uniform distribution of the phases in (B,B + 
A) with the same coupling strength K = 5, K > K c , the 
same m = 6 but different B. 


In the following, physical explanations are given on the 
basis of attraction and repulsion interaction among the 
different oscillators . 

The coupling strength I\ > 0 in Eq. © is attracting 
to synchronization and is repulsive to synchronization 
as K < 0. The incoherent state will stable in the case 
K < 0. Concerning the attractive and repulsive proper¬ 
ties of the coupling parameter, Hong and Strogatz have 
studied the identical oscillators with some couples others 
negatively (the contrarians) and some of positive cou¬ 
pling (the conformist). The contrarians like to be anti¬ 
phase with the mean field and the conformist is easy to 
in-phase. Also other interesting phenomena like travel¬ 
ing wave occurs [H|. There also are references investigat¬ 
ing the phenomena [39|-[44|. In Ref.[39 . the N identical 
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phases with the non-linear coupling is studied and the 
positiveness and negativeness of the coupling parameter 
is controlled by the non-linearity coupling. The dynamics 
of the system is 

Ni 

Bi = (1 - eZf) Y, Kx sin(% - 6 Z ), (10) 

3=1 

where the local order parameter Z\ = Ylf=i e lmej \. 
The repulsive coupling is realized if 1 — eZf < 0. The 
nonlinear coupling will result in phase-locked states, 
while the large nonlinear coupling will give rise to multi¬ 
stable, periodic and chaotic states [39|. In the neu¬ 
ral network, the fast studying model is attributed as 
Oi = Wi + Jf J 2 J =1 Kij sin (Bj ~ 9i) with varying coupling 
strength as I\ij = a cos {0j — Oi). This is actually the 
second Harmonic Kuramoto model. In this way, the at¬ 
tractive and repulsive coupling parameter is achieved de¬ 
pending the difference of the two phases of the two os¬ 
cillators (26j-[i3]. Similarly, the attracting and repulsive 
properties of the coupling strength are the key to explain 
the cluster phenomena in the the higher odder harmonic 
coupling Kuramoto. 

For the parameter K m > K c > 0 , the coupling 
strength can be either attracting or repulsive depend¬ 
ing the difference of the two phases. If all ( 0j — Oi) are 
less than —, then m(0j — Oi) < 2n and they can be 
collectively synchronic and form a cluster and most os¬ 
cillators synchronic in the cluster. Further increase of the 
range of the phases over the oscillators with phases 


greater than — will be repulsed by the oscillators with 
phases less than —, so large phase difference will form for 
these two kinds of oscillators. Because the most oscilla¬ 
tors cluster synchronically, the repulsion to the oscillators 
with phases larger than — is dominant, hence, these os¬ 
cillators could only oscillators with phases large than — 
and the second cluster emerges. All along this way, more 
clusters will appear as the range of the phases becomes 
larger and larger. 

Conclusion and discussion: in generalized Kuramoto 
with the higher order harmonic coupling, the view from 
the symmetry transformation gives the explanation to 
CS both profound mathematical insight and clear phys¬ 
ical understanding. Detailed numerical studies con¬ 
firm the symmetric analysis. The similar analysis could 
extend to the forced Kuramoto model 0 n = ui n + 

YljLi K m smm(0j — 0 n ) + F n (t), with the force tak¬ 
ing the form of F n (t) = F sin fit in the neural learning 
network. Whenever the force is correlated with the oscil¬ 
lators, like F n (t) = F sin(flf — 0 n ), there is no symmetry 
group transformation like Eq. (J4]). Neither can the Ku¬ 
ramoto model with mixed higher harmonic orders cou¬ 
pling have symmetry group transformation like Eq. d4j . 
So new ideas are needed to be explored in the two cases. 
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